WIMP diffusion in the solar system including solar depletion and its effect on Earth 
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Weakly Interacting Massive Particles (WIMPs) can be captured by the Earth, where they even- 
tually sink to the core, annihilate and produce e.g. neutrinos that can be searched for with neutrino 
telescopes. The Earth is believed to capture WIMPs not dominantly from the Milky Way halo 
directly, but instead from a distribution of WIMPs that have diffused around in the solar system 
due to gravitational interactions with the planets in the solar system. Recently, doubts have been 
raised about the lifetime of these WIMP orbits due to solar capture. We here investigate this issue 
by detailed numerical simulations. 

Compared to earlier estimates, we find that the WIMP velocity distribution is significantly sup- 
pressed below about 70 km/s which results in a suppression of the capture rates mainly for heavier 
WIMPs (above ~ 100 GeV). At 1 TeV and above the reduction is almost a factor of 10. 

We apply these results to the case where the WIMP is a supersymmetric neutralino and find that, 
within the Minimal Supersymmetric Standard Model (MSSM), the annihilation rates, and thus the 
neutrino fluxes, are reduced even more than the capture rates. At high masses (above ~ 1 TeV), 
the suppression is almost two orders of magnitude. 

This suppression will make the detection of neutrinos from heavy WIMP annihilations in the 
Earth much harder compared to earlier estimates. 



o 
o 



Of 
i 1 

o 

d ■ 



X 



I. INTRODUCTION 

There is mounting evidence that a major fraction of the 
matter in the Universe is dark. The WMAP Experiment 
gives as a best fit value that Q Ocdm/i 2 = 0.113 ±0.009, 
where ficDM is the relic density of cold dark matter in 
units of the critical density and h is the Hubble parameter 
in units of 100 kms _1 Mpc _1 . One of the main candi- 
dates for the dark matter is a Weakly Interacting Mas- 
sive Particle (WIMP), of which the supersymmetric neu- 
tralino is a favorite candidate. There are many ongoing 
efforts trying to find these dark matter particles, either 
via direct detection or via indirect detection by detecting 
their annihilation products. 

One of the proposed search strategies is to search for a 
flux of high-energy neutrinos from the center of the Earth 
. This idea goes back to Press and Spergel 0] , who cal- 
culated the capture rate of heavy particles by the Sun. 
For the Earth, the idea is that WIMPs can scatter off a 
nucleus in the Earth, lose enough energy to be gravita- 
tionally trapped, eventually sink to the core due to subse- 
quent scatters, annihilate and produce neutrinos. Due to 
purely kincmatical reasons, the capture rate in the Earth 
depends strongly on the mass and the velocity distribu- 
tion of the WIMPs. The heavier the WIMP is, the lower 
the velocity needs to be to facilitate capture. In Ref. 
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[4j, Gould refined the calculations of Press and Spergel 
for the Earth and derived exact formulae for the capture 
rates. In a later paper p|, Gould pointed out that since 
the Earth is in the gravitational potential of the Sun, 
all WIMPs will have gained velocity when they reach 
the Earth and hence capture of heavy WIMPs would be 
very small. However, Gould later realized Q that due to 
gravitational interactions with the other planets (mainly 
Jupiter, Venus and Earth), WIMPs will diffuse in the 
solar system both between different bound orbits, but 
also between unbound and bound orbits. Gould showed 
that the net result of this is that the velocity distribu- 
tion at the Earth will effectively be the same as if the 
Earth was in free space. This approximation is widely 
used today where one further assumes that the halo ve- 
locity distribution is Gaussian (i.e. a Maxwell-Boltzmann 
distribution). 

However, Farinella et al. later made simulations of 
Near Earth Asteroids (NEAs) that had been ejected from 
the asteroid belt. They found that many of these have life 
times of less than two million years. After that time they 
are either thrown into the sun or thrown out of the solar 
system. If this typical lifetime also applies to WIMPs, 
this would significantly reduce the number of WIMPs 
bound in the solar system, as pointed out in Ref. |8|]. 
This in turn would reduce the expected capture and an- 
nihilation rates in the Earth and thus reduce the neutrino 
fluxes. In Ref. ||| , Gould and Alam investigated what the 
implications would be if bound WIMPs would actually 
be thrown into Sun. They investigated two scenarios: 
an ultra conservative scenario where all bound WIMPs 
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are depleted and a conservative scenario where all bound 
WIMPs that do not have Jupiter-crossing orbits are de- 
pleted. In the ultra conservative view solar depletion is 
assumed to be so efficient that no bound WIMPs exist, 
whereas in the conservative view, Jupiter is assumed to 
be faster at diffusing WIMPs into the solar system than 
solar depletion is at throwing them into the Sun. Both of 
these views significantly reduce the neutrino fluxes from 
the Earth for heavier WIMPs. 

However, the truth probably lies somewhere between 
the conservative view and the assumption that solar de- 
pletion is very inefficient, i.e. some WIMPs on bound 
orbits in the inner solar system will survive, but solar 
capture will diminish their numbers somewhat. The aim 
of this paper is to investigate the effects of solar capture 
on the distribution of WIMPs in the solar system and 
the implication this has on expected neutrino fluxes from 
the Earth. We will do this by numerical simulations of 
WIMPs in the solar system and by reanalyzing the pro- 
cess of WIMP diffusion in the solar system. Finally, we 
will apply our results to the case where the WIMP is the 
neutralino, which arises naturally in Minimal Supersym- 
mctric extensions of the Standard Model (MSSM). 

The layout of this paper is as follows. In section [HI 
we will briefly review the history of WIMP capture cal- 
culations for the Earth. In section ITTT1 we will go through 
our assumed halo model and the role of diffusion in more 
detail. In section llVl we will go through the formalism for 
the diffusion caused by one planet and in section we 
add the new ingredient, solar depletion. In section IVT1 we 
present our numerical treatment of the diffusion problem. 
All of this will be put together with the dominant plan- 
ets for diffusion in section rVlII where our main results on 
the velocity distribution at the Earth are presented. In 
the remaining sections we will investigate how this affects 
the capture and annihilation rates in the Earth and will 
present results on the expected neutrino-induced muon 
fluxes in MSSM models in section 11X1 Finally, we will 
conclude in section Ixl 



II. CAPTURE OF WIMPS BY THE EARTH - 
HISTORICAL REMARKS 

Capture of WIMPs by the Sun was first studied by 
Press and Spergel 1985 [3j. Their calculations were 
approximate in nature, especially when applied to the 
Earth. This was refined in a series of papers by Gould 
H H @. In 1987, Gould derived the exact for- 
mulae needed to calculate the capture of WIMPs by a 
spherically symmetric body. When applied to the cap- 
ture by the Sun and the Earth, his approach enhanced 
the capture by factors of 1.5-3 and 10-300 respectively, 
compared to the previous approximations by Press and 
Spergel 0. 

However, in 1988, Gould [5| refined the analysis, tak- 
ing into account that the Earth is well inside the gravita- 
tional potential of the Sun. The velocities of the incoming 



particles are increased when they approach the potential 
of the Sun. This reduces capture substantially. On the 
other hand, bound solar orbits are allowed. Gould real- 
ized that particles scattered by the Earth could become 
bound to the solar system. This scattering can be of two 
kinds: gravitational scattering, which is elastic in that 
the velocity with respect to the Earth is conserved or 
weak scattering off an atom, which can be inelastic and 
either lead to capture by the Earth or make the particle 
bound to the solar system. In this context, an equation 
for estimating the timescales of weak and gravitational 
scattering was developed, following traditions of Opik 

Among other things, Gould concluded that, due to the 
differences in total scattering cross section, the gravity 
of the Earth is more effective in changing the orbits of 
bound particles than is weak scattering. For capture by 
the Earth though, weak scattering is the only process 
at work since gravitational scattering leaves the velocity 
with respect to the Earth unchanged. 

In 1991, Gould continued further, moving his at- 
tention to the gravitational diffusion caused by the other 
planets. Further, he considered the combined diffusion 
effect of Jupiter, Venus and the Earth concluding that it 
will make the velocity distribution isotropic in the frame 
of the Earth. Based essentially on Liouville's theorem, 
this means that the phase space density of unbound and 
bound particles would be the same. Specifically, for the 
most important parts of velocity space, this would hap- 
pen on time scales shorter than the age of the Solar Sys- 
tem. Obviously, such a scenario would substantially en- 
hance capture of heavy WIMPs by the Earth. Further, 
he concluded that weak capture of WIMPs to bound solar 
orbits is negligible, and that one may use the "free-space" 
formulae derived in Ref. for capture, even though the 
Earth is deep within the potential well of the Sun. 

As mentioned in the introduction, the calculations took 
a new unexpected turn in 1999 when Gould and Alam 
interpreted Q the results of Farinella et al. Q- Farinella 
et al. had numerically calculated the fates of about 50 
asteroids of which most were considered to be near Earth 
asteroids (NEAs). They concluded that about a third 
of the considered asteroids will be ejected to hyperbolic 
orbits or, more importantly, driven into the sun in less 
than 2 million years. If the results of Farinella et al. were 
applicable to general Earth crossing orbits of WIMPs, 
the part of velocity space corresponding to bound solar 
orbits would be effectively empty, since the typical time 
scales at which such orbits are populated from unbound 
orbits are generally much longer [6j . The basic results of 
Farinella et al. were later confirmed by Gladman et al. 
[l(j and Migliorini et al. 

To investigate the role of Solar depletion, Gould and 
Alam Q analytically investigated the difference between 
the 1991 case of no solar depletion, and the other ex- 
treme, where there is no dark matter in solar system 
bound orbits at the Earth. In the latter ultra conser- 
vative view, capture is heavily suppressed. For instance, 
they found that WIMPs with masses above about 325 
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GeV could not be captured by the Earth at all. They also 
considered a scenario where WIMPs in Jupiter-crossing 
orbits were not effected by solar capture, the conserva- 
tive view, which they found allows WIMPs up to about 
630 GeV to be captured. (See section ITlI Bl below for a 
more detailed discussion of these cutoff masses.) Both 
of these scenarios significantly suppress the capture rates 
of WIMPs by the Earth and the question to ask is if 
the results of Farinella et al. can really be applied to all 
Earth-crossing WIMP orbits? The orbits of asteroids 
ejected from the asteroid belt are, after all, rather spe- 
cial as they typically arise from resonances. It is thus 
not necessarily so that these results apply to all bound 
WIMPs. We will in the coming sections go through the 
necessary steps to investigate this question in detail. 



distribution pass through the solar system, the velocities 
are boosted and focused by to the gravitational potential. 
At the location of the Earth, the solar system escape 
velocity is V2v e «42 km/s, where we have used the speed 
of the Earth, v 9 ~ 29.8 km/s. Therefore the velocity at 
the location of the Earth, vo, is, according to conservation 
of energy 



w 2 = s 2 + 2vi. 



(4) 



When a spherically symmetric distribution such as 
F a (s) is focused by a Coulomb potential such as that 
of the Sun, the following statement holds: [f| 



F w (w)4irw 2 dw F s (s)4:Trs 2 ds 



(5) 



III. THE GALACTIC HALO MODEL AND 
CUTOFF MASSES 

A. The galactic halo model 

In order to make the calculations concrete, we use the 
Maxwell-Boltzmann model [l^, where the local velocity 
distribution of WIMPs is Gaussian in the inertial frame 
of the Galaxy. At the location of the Sun the distribution 
is 



f v (v)d 3 v 



-v 2 /v 



TV 



3/2,,3 



d s v, 



(1) 







where uo = y |w with v being the three-dimensional ve- 
locity dispersion. We will here use the standard value 
of v — 270 km/s corresponding to vq — 220 km/s. The 
distribution is normalized such that 



f v (v)4:Trv 2 dv = 1. 



(2) 



The velocity distribution can be galileo transformed 
into the frame of the Sun: / s (s), where s = v + vg un , 
and wg un =220 km/s, and averaged over all angles. In 
this special case of a Gaussian distribution the transfor- 
mation can be done in closed formal . As Gould have 
pointed out, the angle between the rotation axis of the 
solar system and that of the galaxy is about 60° which 
makes the velocity distribution very close to spherically 
symmetric, if one considers averages over a galactic year 

200 million years[^. The distribution used is mirror 
symmetric in the galactic plane which means that the 
time of average need only be 100 million years. 

The symbol F s (s) will be used to denote the phase 
space number density 



F s (s) 



(3) 



This can be understood as Liouville's theorem for the 
spherically averaged phase space density, since 



ds 
dw 



F w {w)=F s (s). 



(6) 



Since the velocity w of the halo particles is always at 
least equal to the escape velocity, there will be a hole in 
velocity space so that 



F w (w) = when w < v2t 



(7) 



This is important since capture by the Earth is very sen- 
sitive to F w (w) at low velocities. 

The distribution F w (w) can now be used to calculate 
the distribution as seen from the moving Earth where the 
particle velocity is u = w + v e 



F u (u) =F w (w) =F w (u- v ffl ). 



(8) 



This means that the hole is shifted, so that it is centered 
around — v e . This is visualized by figure^which displays 
a two dimensional slice of the three dimensional velocity 
space. 



B. Cutoff masses when low velocity WIMPs are 
missing 

In the absence of WIMPs gravitationally bound to the 
solar system, the capture by the Earth is totally sup- 
pressed for WIMP masses larger than a critical value. 
To understand this, consider a particle approaching the 
Earth with velocity u at infinity with respect to the grav- 
itational potential of the Earth. If it is to be captured, 
it must be scattered by an atom to a velocity less than 
the escape velocity i>esc & t the atom. By conservation of 
energy and momentum, the particle must have a veloc- 
ity less than (assuming iron to be the heaviest relevant 
element of the Earth) 



where M x is the WIMP mass, and p x is the WIMP mass 
per unit volume in the halo. When the particles of this 



y/M x M Fe 

"cut = 2— in— v esc 



M x - Mi 
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(9) 
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Speed at the Earth (km/s) 

FIG. 1: The ecliptic (cj> — it/2) slice of the particle velocity 
space in the frame of the Earth. The dotted curves show the 
velocity relative to the Earth and the indicated angle 9, is the 
angle of the particle with respect to the direction of Earth's 
motion. The angle <j> determines in which angle we cut the 
velocity sphere. <j) — is the south pole of the solar system 
and <fi = 7r/2 (as shown here) is the slice radially outward from 
the Earth. The region inside the solid semicircle represents 
bound orbits. It's radius is the escape velocity from the Solar 
system at the location of the Earth, but in the frame of the 
Sun. In the same way, the region outside the dash-dotted line 
(an almost perfect semi-circle) corresponds to particles that 
may reach Jupiter. By repeated close encounters with the 
Earth, particles may diffuse along the dotted circles (actually 
spheres) of constant velocity only, keeping u constant, but 
allowing changes in 8 and (j>, as explained in the text. 
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FIG. 2: Cutoff velocity v cnt against WIMP mass M. Only 
combinations of M and « cu t to the left of the line is kinemat- 
ically allowed (in the sense that they can lead to capture by 
the Earth). 



when it approaches the gravitational potential of the 
Earth. Solving for the WIMP mass M x gives 



M 



x,cut 



2wesc («esc + \fu?~ 



u esc) 



(10) 



Here, u is the speed (at infinity) of the approaching par- 
ticle in the frame of the Earth and M x cu ^ is the highest 
allowed mass of the particle if it is to be captured by the 
Earth. The escape velocity varies from 11.2 km/s at the 
surface to 15.0 km/s at the center of the Earth (see sec- 
tion ^^^^ for more information about the Earth model 
we use), and capture is thus easiest at the center where 
the escape velocity is higher. Using v csc — 15.0 km/s, 
we plot in Fig. the relation between u cu t and the cutoff 
mass, M x , cut . 

With Eq. Ijl0|l . we can now relate to the cutoff masses 
in the conservative and ultra conservative views by Gould 
and Alam Q . In the ultra conservative view, we assume 
that only unbound halo particles are captured. Halo par- 
ticles cannot be slower than w cut = 

(V2- 1)« 8 - 12-3 
km/s at, and in the frame of, the Earth (this is also seen 
in Fig. . This gives a cutoff mass of about 410 GeV 
over which capture by the Earth is impossible. This dif- 
fers from the value of 325 GeV for the ultra conservative 
view in Gould and Alam The difference is because 
they used an average escape velocity of 13 km/s instead 
of the maximal one of 15 km/s that we have used in 
Fig. El 

In the conservative view, we assume that Jupiter- 
crossing orbits are filled. This means that all orbits 
outside the dot-dashed curve and the dashed curve in 
Fig. [l]are filled. The lowest velocity WIMP at the Earth 
that is on a Jupiter-crossing orbit is in the lower right- 
hand end of the dot-dashe d curve and it has a velocity of 
""cut = v$(yj2/(l — r & /r % ) — 1) ~ 8.8 km/s (and is mov- 
ing in the same direction as the Earth). r% ~ 5.2r ffl is 
the radius of the Jupiter orbit. This value of u cu t gives a 
cutoff mass of about 712 GeV, whereas Gould and Alam 
p got a cutoff mass of about 630 GeV. The difference is 
again due to the different escape velocities used, but also 
a different velocity to reach Jupiter. We use the value 
it cut ~ 8.8 km/s as indicated above, whereas they used 
an approximation for more general orbits than the one 
giving the cutoff derived here. So, to conclude, in the 
conservative and ultra conservative view, we cannot cap- 
ture WIMPs heavier than about 410 GeV and 712 GeV 
respectively. This is in rough agreement with the results 
of Gould and Alam 0. 

If, on the other hand, the solar system is full of gravi- 
tationally bound dark matter, the velocities can be much 
lower. As the lowest allowed velocity of the WIMPs u cu t 
tends to zero, the mass limit M x cu ^ goes to infinity. 
Typically, most WIMPs in the Galaxy have velocities 
much greater than those of Eq. © , so only a small frac- 
tion of the WIMPs are possible to capture. 

A particle in close encounter with a planet, for instance 
the Earth, may get gravitationally scattered into a new 
direction and a new velocity as seen from the frame of 
the Sun. However, by conservation of energy, the speed u 
with respect to the frame of the planet is unchanged. This 
means that a particle at a particular place in velocity 
space may, by repeated close encounters with the Earth, 
diffuse to any location on the sphere of constant velocity 
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(with respect to the Earth), and nowhere else. 

The location of a particle at such a sphere can be speci- 
fied by the angles at which it passes the Earth. The angle 
is measured between the forward direction of the Earth 
and the velocity vector of the particle, and <j> is the angle 
of rotation around the forward direction of the Earth, 
with cf> — at the North pole of the solar system. 

Figure ^ illustrates how spheres (and circles) of con- 
stant u cross the limit of where particles have bound and 
unbound orbits. This corresponds to the possibility of 
gravitational capture and ejection from the solar system. 

A single planet can diffuse particles along spheres of 
constant velocity only. Therefore, it is clear (from e.g. 
Fig. ^) that orbits with velocities u less than 12.3 km/s 
cannot be populated by the Earth alone. However, since 
the velocity spheres of different planet's are not concen- 
tric (they need not even be spheres when the particles 
reach another planet), the combined effect may diffuse 
particles down to Earth-crossing velocities u less than 
12.3 km/s. This will be investigated in detail in section 
IVIII In order to do so, we must first understand how a 
single planet affects the phase space distribution. 



IV. GRAVITATIONAL DIFFUSION IN THE 
ONE PLANET CASE 

In this section, we investigate the details of what will 
be called gravitational diffusion. We will develop tools 
for detailed investigation of the bound orbit phase space 
density, taking the effects of solar depletion into account. 
We will here start by looking at diffusion effects from a 
single planet only and will take the Earth as an example. 
The exact same formalism is then used for Venus and 
Jupiter as well. 

In this section we assume that when a particle is in 
Earth crossing orbit (perihelion less than the Earth or- 
bit radius i? e and aphelion greater than i? ffi ), long range 
interactions with other planets are less important, and 
can be ignored. This is not a problem, as we in section 
IVIII add the effects of other planets (apart from possible 
resonances). We will in this section closely follow Gould 
0, with some small modifications. 



A. The probability of planet collisions. 

We are interested in calculating the rate at which 
WIMPs with Earth crossing orbits comes into close en- 
counter with the Earth. This will be used to estimate 
how the Earth affects the WIMP distribution. A close 
encounter is an event were the particle's impact param- 
eter is smaller than or equal to some value b max (u). 

Let's imagine the Earth as being spread out on a flat 
ring of inner radius i?, outer radius R + I and thickness 
h, as in figure Now consider a particle with perihelion 
less than the planet orbit radius R and aphelion greater 
than R. Such particles will be said to have Earth crossing 




FIG. 3: The angle of perihelion precession A£, as the orbit 
enters and leaves the disk of the Earth. In this example, the 
plane of the orbit is nearly perpendicular to the ecliptic plane. 



orbits. This is motivated by the fact that due to the pre- 
cession of perihelion, all such orbit ellipses will eventually 
intersect the Earth ring. The small angle the perihelion 
sweeps out, as the orbit ellipse enters and leaves the ring 
is given by 

A£ w tanA£ = 4|tan8i| (11) 
R 

where 8i is the intersection angle between the WIMP 
ellipse and the plane perpendicular to the location vector 
of the Earth R. Since this happens four times during 
each perihelion revolution, the mean probability for such 
a WIMP to intersect the ring of the Earth during each 
WIMP year T x is 

M = W 

The probability for the WIMP to come into close en- 
counter with the Earth is therefore px times the cross 
section a of such an event, divided by the area over which 
the Earth is distributed. However, the length of the path 
which is inside the Earth ring during each encounter is 
oc | cos 02 1 -1 , where 82 is the angle between the axis 
of the ecliptic and u, the velocity of the WIMP as seen 
from Earth [j| . The probability for a reaction with cross- 
section (j, can now be calculated, 

1 a 41 , „ , 1 

'tanOrl — . (13) 



p(a) 



Jcose 2 | 2irRl 2nR 



where we have divided by T 9 to get the probability per 
unit time. The WIMP year can be written in terms of 
u(#, 4>,u) , the velocity of the particle in the frame of 
the Earth, 



T x =(l- 2— cosfl- 



U n-3/2 



T 



u 

and 81 (u) and 82(11) can be expressed in u, 9 and 

R (u + v©) 



(14) 



cos 81 



R ■ v\ 



R\u + v s 



u sm v sm < 



v% + 2uv fS cos8) 1 / 2 



(15) 



n U • VqXR 

cos ©2 = sin 9 cos <p = , (16) 

UVasR 



U 



cot 01 



sin v sm i 



(1 + 2(u/v @ ) cosd + (u 2 /v 2 )(l - sin 2 9 sin 2 (/>)) 



V 



on 



By substituting and rearranging we conclude that the 
yearly reaction probability for an event with cross- 
section a is given by 



p(u, u) 3 cr w 9 



2ttR 2 



7(11) , where 



(18) 



7(u) 



3 7rsin 2 | sin0cos0| (1 — 2u/v fB cos8 — u 2 /v 2 ) 3 / 2 n \ 
% (l + 2(u/v m )cosd + (u 2 /v 2 )(l-sm 2 " ' " ' ^ 2 " 



'sin 2 0)) J 



Eq. Ijl9(l was first derived by Gould [S| in a very similar 
way. Among other things, he used it to calculate the 
"typical timescales" at which particles diffuse between 
different velocity space regions in the absence of solar 
depletion. It is also used for calculating the probability 
of weak scattering of WIMPs at the Earth. 

The equations above are derived under some (geomet- 
rical) approximations with the aim of getting the correct 
scattering probabilities on average. There are however a 
few pathological cases where the geometrical model used 
above breaks down. This happens when </> = 0, </> = 7r/2, 
8 = and 8 = n, in which case the probabilities above 
are unphysical. Since this only happens for these few 
special cases we will artificially solve this by adding a 
small angle (of about 1 degree) to 8 and <j> when close to 
these regions. Note that in principal, the problems could 
be resolved, by making A£ a function of the full set of 
orbit parameters, but this is unnecessarily complicated 
for our purposes. For the interested reader, we refer to 
a detailed investigation of the mathematical properties 
of 7 as presented in and @ . To test our solution of 
adding a small angle in these pathological cases, we have 
investigated the effect of further increasing the small an- 
gle added and conclude that the actual value chosen is 
not important for the final results. This is reasonable, 
since orbits in the vicinity of these critical regions are 
quickly deflected into other orbits anyway. 



B. Gravitational scattering on a planet 

Now that we have learnt how to calculate the proba- 
bility for particles to come into close encounter with a 
given planet, it is time to apply this to gravitational dif- 
fusion. For the Earth, we were mainly interested in those 
particles crossing the sphere of one AU during each revo- 
lution, since they have a chance of hitting the Earth (and 
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FIG. 4: Scattering off the Earth in velocity space. A fixed 
impact parameter fixes 5, but r\ is evenly distributed. 

possibly be weakly captured to it) within each perihelion 
precession revolution. 

The gravitational scattering probability is dependent 
on the angular distance between the velocities before and 
after scattering: u and u', such that small deflections are 
more common. The angle can be related to the impact 
parameter 6, 



5(b) 



2 arctan 



bu 2 
MG' 



as well as 



5(u', u) = arccos(u' • u). 



(20) 



(21) 



where ~ denotes unit vectors. The scattering angle above 
is the one given by Rutherford scattering (see e.g. [l^)- 
Gould used an approximate formula when deriving the 
typical time scales 5(b) = R^v 2 sc /(bu 2 ). The two 
differ at very small impact parameters, and we use the 
full expression in our calculations. 

As mentioned before, scattering can only change the 
direction and not the velocity, and we are therefore deal- 
ing with random walk on spheres of constant u. The 
direction r\ of the scattering is evenly distributed, as seen 
in Fig. ^ where the scattering setup is shown. The arc 
length is fixed by 5(b), but the scattering direction is 
evenly distributed. 

The cross-section for scattering between 5 and 5 + dS 
is da = 2Trbdb, so the yearly probability for scattering in 
this range is, (using Eci. lTB|) 



dp(u, b) 3 2nb % 
~dbT~ ~ 2,'ttE 2 ~u, 



T(u)- 1 . 



(22) 



This can be rewritten in terms of the scattering angle 
5(xx',u), 



d P (u,b(S)) = 

dsn 

3 2n Vq 



2ttR 2 



7 (u) 



_i (Mj,G) 5 



(5(u',u)/2) 



2sin d (£(u / ,u)/2) 



(23) 
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FIG. 5: An example of the timescales of particle scattering. 
The color bar indicates the time for which there is a 10% 
probability of scattering an angle S = n/2 ± 7r/64, depending 
on the present location of the particle. In general, time scales 
are shorter at lower velocities. By repeated close encounters 
with the Earth, particles may diffuse along the dotted circles 
(actually spheres) of constant velocity only, keeping u con- 
stant, but allowing for changes in 6 and (j>. The figure shows 
the cj> = 75° slice of the particle velocity space only. 
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The right-hand side of Eq. 1)23(1 may look like a negative 
probability density, but this is artificial since integration 
should be done for decreasing S's. We integrate Eq. lt2*5|) 
analytically and use that expression whenever numerical 
values of the scattering probability to <5 + AS are needed. 

To get a feeling for the significance of the diffusion, 
we solve Eq. <|23[1 to obtain the typical time scales for 
scattering a given angle to occur. As an example, we look 
at the time scales for which the probability of scattering 
with S = tt/2 ±7r/64 is 10%. This is illustrated in Fig. El 



C. The bound orbit density and orbit capture from 
the halo 

Let us now define the bound orbit density, n(u) to be 
the number of bound particle orbits per infinitesimal ve- 
locity and solid angle on the velocity sphere. The orbit 
density is thus free from information about the particle 
location along its elliptical orbit. The total number of 
bound particle orbits in a thin shell of radius u is 



dN = duu 1 



Sl=bound orbits 



dVt n(u) 



(25) 



We will now divide each velocity sphere into cells (that 
can at this point be viewed as infinitcsimally small). The 
number of particles scattered between two locations on 



a sphere of constant velocity in a given time must be an 
integral over the source cell i and the destination cell j, 



dN r 
~dT 



dP({3, a) 
db 
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In our case, the destination space is conveniently spanned 
by the scattering angles 8 and r\. The density of bound 
orbits scattered from cell i to cell j evolves with time as 



driji 
~dT 



dil 



dn dp(u, S) 

K dS 27 ^hTT" {xu - (26) 



where Kj is defined to be the region in S-r)-sj>&ce cor- 
responding to scattering from the i to the j cell. The 
scattering probability to the [S, 8 + AS] band is evenly 
distributed over all cells in that region. Numerically this 
is implemented as as loop over 5 as measured from the 
center of the source cell. The probability is then dis- 
tributed over all discrete cells whose centers are inside 
the current band. 

It is important to understand that what we are consid- 
ering is the movement of the particle orbits, as opposed 
to the particles themselves. This means that we do not 
need to calculate the actual particle trajectories. When 
we are interested in the actual particle densities, we pick 
orbits from the orbit densities. Note that there are two 
points on the orbit that pass a given radius, but due to 
the perihelion precession, any given particle could show 
up at any of these orbit locations depending on the angle 
of perihelion. Since we have anticipated mirror symmetry 
in the plane of the solar system, each particle is smeared 
out at four indistinguishable orbit locations on the sphere 
of constant u. This can always be done, regardless of the 
existence of symmetries in the free distribution. If the 
free distribution is not mirror symmetric in the ecliptic 
plane, it can be forced to have this (in this case artificial) 
symmetry by averaging, as long as we are only interested 
in the absolute capture of WIMPs in the Sun or planets. 

The equations derived above apply only to particles 
which are already gravitationally bound to the solar sys- 
tem. We now turn to the calculation of the bound orbit 
density capture rate] Anjf/T m from the distribution of 
free particles. We will use the local phase space density 
Ff(u) [Particles/ (m 3 m/s)]. 

Consider the distribution of particles Ft (u) passing the 
Earth with impact parameter b. The number of particles 
scattered an angle 8{b ± db/2) in a given period of time 
T, is 



Tu2nbdb F f (u). 

dV 



(27) 



According to Eq. I|2U|I they are scattered an angle S(b). 
Using the relations (|24|l we conclude that the bound orbit 
density at the cell j will evolve with time as 
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caused by gravitational scattering from the halo. We 
now have equations for gravitational diffusion as well as 
capture to the solar system. 



D. Relating the phase space density F(u) and the 
bound orbit density n(u) 

The ideas of the last section can be used to write down 
an expression for the phase space density, which is what 
we need for the weak capture calculations. The relation 
between the phase space density F(u) and the bound or- 
bit density; n(u) is derived as follows. 

For a given orbit in the population of bound orbits, 
we use Eq. Ijl8|) to calculate the number of orbits that 
will pass trough an area a each year. We now consider a 
volume dV in space with base area a and height h such 
that h is parallel to u. A particle passing trough the area 
will spend a time h/u in the volume. This means that 
the fraction of the WIMP year spent in the volume in 
case of an encounter is 



uT x 

The fraction of orbits passing through a each WIMP year 
is 

T e X 

Therefore, since F (u) is the number of particles per 
du 3 dV, the relation between -F(u) and n(u) is 



F{u)dV = n{u) J^P±l^ Tv 



(29) 



or 



F(u)dV = n(u) h 



P(°-, u) 



3 dV v 9 -l 
2 -kHt u z 



One should note that by construction F(u) above is valid 
in the frame of the planet. However, the right hand side 
of the equations above presumes the planet to have a 
constant velocity during the encounter so that du 3 are 
equal to the velocity volume element in the frame of the 
Sun, du> 3 . 



V. SOLAR DEPLETION OF BOUND ORBITS 

In the previous section, we investigated the evolution of 
the bound orbit densities due to scatterings from other 



bound orbits and from free orbits. We have one main 
piece remaining to be studied, and that is the effects of 
solar depletion, i.e. how much of the bound WIMPs are 
actually captured by the Sun, thus reducing their density 
in the solar system. 

We have done this by numerically calculating the ac- 
tual motion for different WIMP orbits in the solar system 
over 49 million years. As a measure of the quality of the 
numerical methods, we have also calculated the fates of 
the 47 asteroids studied by Farinella et al. 0, as pre- 
sented in appendix lAl 



A. The numerical methods and conditions 

We have numerically integrated the orbits of about 
2000 particles in typical Earth crossing orbits in order to 
estimate the solar depletion. The particles were spread 
out on the bound velocity space with random initial po- 
sitions on the Earth's orbit. We have mainly used the 
Mercury package 0] by Chambers for the integration. 
It has the most important numerical algorithms, such 
as Everhart's 15:th order Radao with Gauss-Radao 
spacings, and the equally well-known Bulirsch-Stoer 0] 
algorithm. Both are variable step size algorithms dedi- 
cated to many body problems, and are commonly used in 
asteroid research for problems similar to ours. The pack- 
age also includes a set of symplectic algorithms, which 
have been used for some tests. By looking at some test 
orbits, we found that the symplectic algorithms (at least 
as implemented in the Mercury package) were slower and 
less accurate for our setup. The tested symplectic al- 
gorithms were "MVS: mixed- variable symplectic" 0] as 
well as "Hybrid symplectic/Bulirsch-Stoer" [l4| . 

The calculations included the test particles, the Sun, 
the Earth, Jupiter and Venus. Other planets were not 
included as they are believed to be sub-dominant. The 
Bulirsch-Stoer algorithm was used to calculate the orbits 
of all test particles, as well as the planets, during a time 
of 49 million years. This took about 35 000 CPU hours, 
on a variety of Linux and Alpha machines. A wide range 
of different accuracy parameters were used, from 10~ 14 to 
10~ 8 , to evaluate the role this plays. The numerical rep- 
resentation of the real numbers limits the benefit of going 
past about 10~ 12 . The final choice of 10~ 10 is a balance 
between time and accuracy. A recent publication |18| in 
the subject of numerical simulations of a special set of 
Jupiter crossing asteroids, came to a similar conclusion; 
When using the Bulirsch-Stoer algorithm for their calcu- 
lations, they found accuracy parameters in the range of 
10 -9 -10 -8 to give statistically similar results as 10 -12 . 
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In the comparisons carried out, this gave results very 
similar to those with higher accuracy parameters. The 
comparison with the Radao algorithm gave qualitatively 
similar, however not identical, results with a similar cal- 
culation speed. In some occasions however, the Radao 
algorithm gave a higher solar depletion for particles with 
very high velocity (relative to the Earth), u > 50 km/s. 
This is not of much concern for our purposes though as we 
are mainly interested in much lower velocities for Earth 
capture to be efficient. 

For ordinary asteroid calculations, a point mass ap- 
proximation combined with collision detection is suffi- 
cient. Our case is a little more delicate since WIMPs 
may pass through the planets. To handle this, the grav- 
itational routines were modified to use the real gravita- 
tional potentials inside the planets. 

For Jup iter and Earth, we used "true" mass distribu- 
tions 19, 201. For Venus we rescaled the mass distri- 
bution of the Earth and removed the liquid iron core. 
Other planets included in tests where assumed to be ho- 
mogeneous. The improvement allows the particles to pass 
through the planets without being infinitely scattered by 
a point mass, making the calculations more realistic and 
numerically stable. For completeness, it would be inter- 
esting to add more planets to the simulations, but it is 
unfeasible to do as it slows down the calculations too 
much. We also believe, that we have included the most 
important planets in our simulations. 



B. The results of the numerical simulations 

The solar depletion was mainly calculated for par- 
ticles in eight planes of u space, with the <j> values 
0,15,30,45,60,75,90 and -30 degrees (the = -30° 
plane was used to investigate the expected radial mirror 
symmetry of the results). Our solar depletion results are 
not as bad as Gould feared || ; Most of the particles sur- 
vived two million years. Nevertheless, solar capture is 
too large to be ignored. Figs El and show the <j> = 75° 
plane, and the times after which the particles hit the Sun. 
We note that ejection is much more common at Jupiter- 
crossing orbits. This is in compliance with the fact that, 
according to the scattering model used here, the proba- 
bility of scattering for such orbits is high. The fact that 
there is a large region at —50 km/s where there are no 
ejections or sun captures, is in agreement with the qual- 
itative results by Gould, presented in his 1991 paper Q 
(see his Fig. 3, where he assumes that the filling times 
are about the same as the time of ejection). Apart from 
the calculations shown here, some extra calculations were 
carried out for relative velocities lower than 15 km/s. The 
results of those calculations were incorporated and used 
in the same way as the others. 

Another important, however simple, result is that 
there seem to exist a mirror symmetry in the in-out direc- 
tions. This is expected, since particles may hit the Earth 
both on its way out and on the way back on its perihe- 
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FIG. 6: The time (linear scale) for ejection (blue/dark gray) 
and capture in the Sun (yellow/light gray) of a set of test par- 
ticles. Each bin represents only one particle, so the statistical 
error is high. However, this figure is typical for all angles, 
except that the plateau of fast solar depletion at large "back- 
ward" velocities are raised when <j> approaches 90°. Some 
particles survived in the Solar system for the whole of the 
simulation. Those particles are marked with black dots. 



lion revolution. Considering Fig. [3 it's evident that this 
is equivalent to a symmetry in the sign of <p claimed by 
Gould; that the <f> and — cases are identical. 

The particle orbits were evenly distributed in velocity 
space, but we solve the diffusion equations on spheres of 
constant m, hence we interpolate our results. What we 
need to extract from our numerical simulations is the de- 
pletion frequency, i.e. the expected depletion probability 
per given time. Since the form of the actual distribu- 
tion, of which the results of the numerical calculations 
are samples, are unknown, the most reasonable way to 
estimate the depletion probability per unit time is 



fsv 



1 



Tsu 



(31) 



Figure [S] shows the logarithm of 1/ f, interpolated onto a 
sphere of constant u, namely u = 40 km/s. 



VI. THE EVOLUTION EQUATIONS FOR ONE 
PLANET 

In the previous sections we have presented the analytic 
expressions for the scattering of bound orbits to other 
bound orbits, Eq. J2EJl, as well as capture from free to 
bound orbits, Eq. I)28[l. We have also, by numerical sim- 
ulations, estimated the rate at which orbits are sent into 
the Sun and thus captured. We are primarily interested 
in how the bound orbit density evolves with time, and 
will here write down the dynamic equations in a form 
suitable for numerical work. 
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FIG. 7: The time (log scale) for ejection and capture in the 
Sun of a set of test particles. This figure is identical to figure 
El except for the time scale which is logarithmic. In this 
scale, it's easier to see that there is a small region at —30 
km/s where the solar depletion occurs directly. This is not 
surprising, since this region corresponds to particles with very 
low velocity in the frame of the Sun. The plateau of direct 
solar capture extends further in the special case of (j> = 90 
(not shown) which allows extremely elliptic, or radial orbits. 
(The plane of start positions is then parallel to the ecliptic 
plane.) 
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FIG. 8: The solar depletion at the it = 40 km/s sphere. The 
color bar indicates the logarithm of the typical depletion time 
l//sun- The region to the right are the free orbits, for which 
the solar depletion is irrelevant. At a "backward" velocity of 
30km/s, the Sun-depletions is greater, in agreement with the 
previous figures of this section. In understanding this figure, 
it may help to take a look at the u — 40 km/s line of figure^ 
which corresponds to the central horizontal (<j> = 90°) plane 
of this figure. 



The bound orbit density develops in time in the fol- 
lowing way, 



au u = au u 
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, (32) 



where rij is the number of orbits in the small cell |3l| 
j of the sphere. The sum over i is the flow from and 
to the other bound cells. The rijf and n S j terms are 
representing capture from unbound orbits and capture of 
bound orbits by the Sun, while the nfj term represents 
the ejection of bound particles. 

We will now reformulate Eq. I|32[l in matrix form suit- 
able for numerical calculations. Let us first define our 
state vectors, 



X 




(33) 



where N s is the number of particles captured by the Sun, 
ni is the bound orbit density and Ff is the velocity num- 
ber density of free (unbound) orbits. If the cells i are 
small enough, the various densities can be considered 
constant over each cell. Using this and the fact that 
the n part of the integration is independent of n(u) and 
Ff(u), this means that Eqs. l|26|l and (|28|l can be written 
as 



dnji 
dt 

dn 3f 

dt 



p°*°ni{u) and 



The solar capture can be written in the same way: 



dn s { 
~dt 



: P s >i(ii) 



(34) 
(35) 

(36) 



The p:s can be considered as the probability per unit 
time to transfer particles/orbits from and to the various 
cells. A positive p means that we transfer to the cell 
and a negative p that we transfer from the cell. The p°^° 
element requires an explanation. This is the probabil- 
ity per unit time that an orbit in cell i is not scattered 
to another bound or free cell, i.e. this term includes all 
the scattering out to both other bound orbits, and un- 
bound orbits. The probability for solar capture though is 
handled separately by p^. As the various entries in the 
state vector X have different units (N s is a number, 7ij is 
the the orbit density and Ft is the number density), the 
required conversion factors are also included in the p:s. 
The cells can be of various size, and these sizes are also 
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included as weights in the p:s. We will not write down 
explicitly the expressions for the p:s as they are found 
elsewhere: p bb and p b ' can be extracted from Eqs. I|26(l 
and (|28[) respectively, while p sc , on the other hand, we 
extract from our numerical simulations of solar capture, 
discussed in section Ivl E.g. the p sc :s for the cells on the 
40 km/s sphere can be read off from Fig. [S] 

Integrating Eqs. I|34|) - (|36[l over time, and replacing dt 
with a discrete time step Ai, we can write the evolution 
of the state vector X as 

X(t + At) = T(At)X(t ) , with 

/l P sc \ 

T(At) =01 + P bb (l - P sc ) - P sc P bf (37) 

\0 1 J 

The first row describes the solar capture of bound orbits. 
The second row describes the development of bound or- 
bits, and the capture of free orbits. Its height is given 
by the number of bound cells. The last row is a little bit 
special. One may propose that bound WIMPs scattered 
to unbound orbits should give a contribution in the sec- 
ond column. However, such particles will not meet the 
Earth again, so the lowest part of the matrix should only 
do the job of keeping the unbound phase space density 
constant. The size of the last row unit matrix is of course 
given by the number of free orbit cells Ff. The matri- 
ces P can be regarded as transition probabilities (for the 
given time interval At). The elements in the P matrices 
are given by 

P%(At) = pf f At (38) 
P^(At) = pfAt (39) 
Pff(At) = d^At (40) 

for free to bound orbits, bound to bound orbits and solar 
capture respectively. 

B. The bound orbit density at arbitrary times 

Eq. I|37|l describes the evolution of the state vector X 
during a time step At. We can write the time develop- 
ment operator that takes us to any time t as 

U(t) = [T(Ai)] t/At X(t +t)= V(t)X(t ) (41) 

The exponentiation of T can be done either by diag- 
onalizing T, or (for applicable times t) by repeatingly 
quadrating T. We have calculated and diagonalized T's 
with a variety of different cell configurations. A simple 
polar grid is a good first choice, but it has a large spread 
in shape and area of the cells, which means that valuable 
memory and calculation time is wasted. Therefore, we 
have used cells with the shape of spherical triangles, built 
from icosahedrons or octahedrons. The cells of the body 
were successively divided in four nearly identical spheri- 
cal triangles, until the right number of cells were reached. 




FIG. 9: The discretization of space were made using trian- 
gles. The number of cells on the displayed triangles are 512, 
1280 and 2048, of which the last one was used in our final 
calculations. 



The velocity space of each planet was built up of about 
65 spheres, usually with 2048 cells each, which means a 
total of about 130 000 discrete cells for each planet. 

If the octahedron is used as a starting object, it's pos- 
sible to rotate the sphere to obtain mirror symmetry in 
the in-out (radial in the solar system) and up-down di- 
rections. Since the problems to solve possess the same 
symmetries, this reduces the size of the state vectors by 
a factor of four, and the time evolution operators by a 
factor of 16. 

Most of the numerical calculations take place in this 
compressed space. By making this a run-time option, 
we have verified that this does not introduce any errors. 
Great efforts have been put in verifying the consistency 
of the time evolution operator. As a simple example, the 
probability for a particle to end up anywhere is unity. 
Making use of the mirror symmetry of the equations, it 
has been possible to calculate and diagonalize P's with up 
to some 10 8 elements. The set of time evolution operators 
of a planet, can be thought of as one huge block diagonal 
time evolution operator for the 130000-dimensional cell 
space. The block shape of the matrix is dependent of 
the conservation of energy when a particle is repeatingly 
scattered by a single planet. We have also checked the 
robustness of our results as the number of cells varies. 
If one uses too few cells, one would expect that the ef- 
fect of diffusion is underestimated as small-angle deflec- 
tions (smaller than the cell size) are then artificially sup- 
pressed. At velocities above 8 km/s, our results do not 
change significantly when going from 1024 to 2048 cells. 
Below 8 km/s, however, the resulting WIMP density is 
somewhat larger in our simulation with 2048 cells than in 
our simulation with 1024 cells. This would indicate that 
the density at these low velocities could go up somewhat 
if we used even more cells. However, it is not possible 
to increase the number of cells further, as it is only fea- 
sible to perform these simulations when the full velocity 
space can be maintained in the computer memory simul- 
taneously. It is also not reasonable to perform this part 
of the calculation more accurately than other parts, like 
the solar capture discussed in the previous section. 

We have now set up a framework for diffusion from 
one planet. We have done this following the scheme set 
up by Gould HJ, with some small modifications and im- 
provements. Our main goal has been to make it possible 
to include the effects of solar depletion, and hence we 
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have formulated the diffusion problem in a form suitable 
for numerical work, where the inclusion of solar capture 
is easily done. 

In the next section we will put all of this together, 
where we also include the diffusion effects of the other 
(dominant) planets. 



VII. THE VELOCITY DISTRIBUTION AT THE 
EARTH: COMBINING THE EFFECTS OF 
JUPITER, VENUS AND THE EARTH 
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We have so far considered the diffusion caused by one 
planet at a time and the effect of solar capture. We are 
now ready to include more than one planet in our treat- 
ment. In section Hvl where we investigated the diffusion 
effects caused by one planet, we saw that one planet can 
only change the direction and not the velocity of a WIMP. 
However, WIMPs that have different directions, but the 
same velocity at one planet, will not only have different 
directions, but also different velocities at another planet. 
Hence, the main effect of including more planets in the 
diffusion is to diffuse particles also to different velocities. 
We thus have a mechanism to populate a larger part 
of the phase space at Earth, and this process is hence 
very important, especially for heavier WIMPs. We will 
here include the diffusion effects of Venus, the Earth and 
Jupiter as these are the planets dominating the diffusion 
mechanism 6l. 



Transformation of coordinates and bound orbit 
density when changing planet 



FIG. 10: A detailed view of the <f) = 75° slice of the particle 
velocity space at the frame of the Earth. The solid archs 
represent particles in Venus-crossing orbits, and their speed 
at the location and in the frame of Venus. Since the archs 
are part of spheres of constant velocity at Venus, they show 
the possible directions of diffusion caused by Venus. Jupiter 
crossing orbits (not shown) fill the region between the dash- 
dotted and the solid semicircles. The low velocity region of 
this figure will be discussed around figuretTTlin section IVlI El 



The change of frame consists of two Galileo transfor- 
mations, as well as the journey of the particles in the 
potential force of the Sun. Since the first is just a change 
of origin in the 6-dimensional phase space, the change of 
frame obeys Liouville's theorem, 



F ? (u s ) =F e (u e ) 



(48) 



Using Eq. (j20J, the orbit densities at the two locations 
can now be related as 



The velocity and angles in a planet-based coordinate 
system at a planet with orbit radius a and velocity v, can 
be converted to the coordinates of another planet via the 
energy, angular momentum and inclination. This is not 
enough for the specification of the exact location of the 
particles, but we are only interested in the shape and 
orientation of the orbits, 
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As an example, we will transform the various densities, 
as seen in the frame of the Earth to the corresponding 
quantities at Venus. 



n s (u s ) = n ffi (u e (u s )) Jtst (u 9 ), where 
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(49) 



(50) 



is the Jacobian. 

Using these transformations, it is possible to investi- 
gate how a sphere of constant velocity at a specific planet 
will look when the particles pass the Earth. Figure ITU1 is 
an example of this. 

The archs of constant velocity in the frame of Venus 
are shown to indicate the directions of diffusion caused 
by that planet. Since the lines of constant velocity at 
Venus cross the u e =12.3 km/s line, Venus may diffuse 
particles into the important u ffl <12.3 km/s region. 



B. Solving the many body diffusion problem 

Consider a point in velocity space, in the frame of the 
Earth, and in Jupiter- crossing orbit. At this point, the 
bound orbit density n(u) takes on a value ua at a given 
time tg. Call this density, transformed into a specific 
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point in the reference frame of Jupiter ub- Now, af- 
ter a short period of time, from now on called step size, 
the Earth may have increased (or decreased) ua by an 
amount driA, and Jupiter may have increased (or de- 
creased) ub by dris- Since ua and ub are really a mea- 
surement of the same density, and there are two processes 
(interactions with the two planets) affecting the the dif- 
ferences, the orbit densities after the time step are given 
by 

ua — > n' a = tia + driA + JdriB 

(51) 

ub — > n' B = us + dris + J driA 

where J is the Jacobian for the transformation. Note 
that the step size introduced above is the step size after 
which transfer of densities between planets occur. For 
the diffusion effects of the individual planets during this 
step size, we use much smaller time steps. 

In order to transfer the orbit densities from one planet 
to another in a numerically reasonable way, all cells at 
each planet is matched to the correct cells on the other 
planet. Since there is not a one to one correspondence 
between the cells of different planets, we need to inter- 
polate between cells. We use a linear interpolation, but 
have also checked that a simpler nearest neighbor inter- 
polation gives similar (but more noisy) results. 

The velocity spaces of all pairs of involved planets were 
tessellated, in order to create the matrix of linear inter- 
polation. This means that each cell was identified to 
constitute the corners of to up to six octahedrons. All 
transformed points were identified to belong to a single 
octahedron, and the location of the transformed point 
was given as a linear combination of the octahedron cor- 
ners. This linear combination was then used as interpo- 
lation for the densities. 



C. Numerical issues 

In the previous subsection, our scheme for taking care 
of the diffusion effects of more than one planet was out- 
lined. We will here discuss the measures we have taken 
to make sure that our numerical implementation is stable 
and does not introduce numerical artifacts. 

In order to further improve the stability of the inter- 
polation between the planet cells, the orbit densities n 
are never interpolated directly. Instead, all interpola- 
tions are done between phase space densities F, and then 
converted to the n-space of the respective planets. The 
phase space density F is a slowly changing function, while 
n is not. This is so since among other things, the rough- 
ness of 7, Eq. 1)19)) . is inherited to n, but removed again 
when F is calculated. 

At any time, the densities at the two planets must be 
consistent with each other so that a density at a partic- 
ular point in one frame matches that of the point trans- 
formed to the other planet, as described by Eq. (|49|l . 
Small interpolation errors can build up with time though, 



and we need to take care of this potential problem. To 
force the densities at the two planets, ua and tlb to be 
consistent, they were regularly averaged as follows: 

(n' A + Jn' B )/2^ n' A 
{n' B + J- 1 n' A )/2^ n' B 

From an analytical point of view, this is not needed, but 
it turns out to be a good way of making the algorithm 
more numerically robust. 

The results are stable with respect to step size as well 
as shape of the velocity space used. This is particularly 
true when the averaging outlined above is done. It is not 
necessary to perform the averaging after every time step, 
instead we can perform it much more seldom. Even if we 
perform the averaging for all velocity spheres, it turns out 
that it is unimportant in the region above u «10 km/s, 
where the processes are slower and stable anyway. In 
the steep region below u — 7 km/s, averaging is needed 
though to keep the stability. We have verified that in 
the limit of very small step sizes, the unaveraged results 
approach the ones with averaging even in this region, but 
averaging allows us to get better accuracy and stability 
even with longer step sizes. We have also verified that 
the results are quite stable with respect to the averaging 
frequency. The result figures of the previous section show 
the results of a small step size; 16 thousand years. 

A related problem is that even though the Jacobian 
determinant of Eq. (|50l) is mathematically valid, linear 
interpolations do not assure conservation of mass. This 
means that when repeatingly transferring density infor- 
mation between a pair of planets, one can not be sure 
that the interpolation does not, in error, introduce or re- 
move mass from the system. These artificial 'sources' or 
'sinks' need to be removed. While it is not possible to 
do this on a cell by cell basis, we have investigated the 
total mass transferred to and from each planet and used 
this to renormalize the mass transferred to ensure mass 
conservation. In equilibrium, the error is quite small; un- 
der one percent, but when a distribution is built up, the 
error can be larger than that. 

We have further tested that the step size is not critical 
for the results. This indicates that our numerical imple- 
mentation, with the stability measures outlined above, is 
stable and that the possible errors are under control. 

We have now presented two methods to keep the possi- 
ble numerical artifacts under control: the averaging pro- 
cess to keep the densities consistent between the planets 
and the renormalization process to force mass conserva- 
tion. Both of these make our algorithms both more stable 
and reliable. 



D. Investigation of Jupiter— crossing orbits 

It turns out, that the density of Jupiter-crossing orbits 
is independent of the diffusion effects of the Earth as well 
as those of Venus. This is expected, since the mass of 
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FIG. 11: The </> = 75° slice of the particle velocity space at 
the frame of the Earth, the low-velocity region. The region 
outside the solid u ~ 9 km/s line is populated by the combined 
effect of Jupiter and the Earth. The diffusion effect of Venus 
is needed for particles to reach the region between the 9 km/s 
and 2.5 km/s lines. 



Jupiter is so much larger, and the scattering probability 
increases with the planet mass squared, see Eq. (|23|) . 

To investigate this, we have numerically solved the 
Earth-Jupiter diffusion system in two ways: calculat- 
ing the evolution with Jupiter alone, as well as solving 
the two body diffusion problem with the methods de- 
scribed above. In either case, it takes only a couple of 
million years for Jupiter's Earth crossing orbits to come 
into equilibrium with the unbound orbits. This means 
that for Jupiter-crossing orbits we can safely neglect the 
diffusion effects of the other planets and let Jupiter fill 
these orbits alone. It also turns out that the diffusion of 
Jupiter-crossing orbits is so much faster than solar de- 
pletion, and we can thus ignore solar depletion for these 
kind of orbits. 

We can then already now see that the ultra conserva- 
tive view in Gould and Alam is too pessimistic and 
that at least as many bound WIMPs as in the conserva- 
tive view remains in the solar system. We will next see 
what the fate is for bound orbits further inside the solar 
system. 



E. Investigation of the Earth— Venus— Jupiter 
system 

Inspired by the last subsection, we will from now on 
keep the density of Jupiter crossing orbits constant and 
focus on the combined diffusion effects caused by Venus 
and the Earth. The locking of the Jupiter crossing orbits 
is done in the same way as for the free orbits, see Eq. (|37|) . 
with the forced insertion of an identity matrix in the time 
evolution operator. As mentioned above, this is justified 
by the fact that diffusion of Jupiter-crossing orbits is so 
fast that we can view these orbits as constantly being 
filled from the halo. We also change the interpolations 
between the Earth and Venus so that Jupiter-crossing 
orbits are excluded (as they are filled by Jupiter). 



Before going through the results, let us spend some 
time going through the diffusion processes in the low ve- 
locity region (as seen from the Earth). Figure ITT1 is a 
closer view of the space of low-velocity orbits. If we ig- 
nore the filling effects of Jupiter, the Earth would have 
to diffuse WIMPs all the way from the unbound orbits, 
starting at the u = 12.3 km/s sphere. They could even- 
tually reach the Venus-crossing orbits to the left of the 
figure. Venus could then act to diffuse the particles along 
its spheres of constant velocity. It is evident from the 
figure that the combined effect of the Earth and Venus 
could possibly populate all orbits outside the u = 2.5 
km/s line. By numerical simulation of the Earth- Venus 
system alone, it turns out that solar depletion is so strong 
that gravitational diffusion can only make a small con- 
tribution to the particle density below 12.3 km/s. This 
is no big surprise, since comparing the scattering times 
in Fig. [5] and the solar depletion times of Fig. [7| we see 
that solar depletion is indeed very strong. 

If we instead use the knowledge about the density of 
Jupiter crossing orbits, the situation is very different. 
The Earth can scatter particles directly from the bound 
Jupiter crossing orbits, starting at u ss 8.8 km/s, as op- 
posed to 12.3 km/s for free orbits. Furthermore the time 
scale of scattering, as well as the angular path the WIMPs 
have to travel is much shorter, especially in the low veloc- 
ity region. Hence, solar depletion will not be as effective 
when we include Jupiter, as the time scales for diffusion 
are more comparable to the solar depletion time scales. 

In our full calculations, we will (as mentioned above) 
keep Jupiter-crossing orbits fixed and include Venus and 
the Earth in the diffusion process. The calculations start 
with a Solar system empty of dark matter, five billion 
years ago. The step size (that is, how often the diffusion 
effects arc added to the other planet) was usually some 
hundred thousand years. The first ten million years were 
typically calculated using smaller step sizes, such as 10 
thousand years. The densities converge to their final val- 
ues within a time of 500 million years. An example of 
the resulting phase space density at a sphere of constant 
velocity is given in Fig. 1121 It is important to remember 
that the free distribution was averaged over a period of 
100 million years. After such a time, the bound densi- 
ties take on their final values within about 25%, which 
is an indication that the results might vary slightly dur- 
ing the galactic (half) year. In practice, this has little 
effect, since the typical time scales for equilibrium (see 
section IV 111 Bll between capture and annihilation in the 
Earth are much longer than that and will average out 
these small variations over the galactic (half) year. 

The resulting velocity distributions for the slowly mov- 
ing particles are shown in Fig. 1131 The ultra conservative 
and conservative curves represent the contributions from 
unbound, as well as unbound plus Jupiter-crossing orbits 
respectively. For these, we again see the cutoff velocities 
of 12.3 km/s and 8.8 km/s as explained in section UlI Bl 
The result of our full simulation, but ignoring solar de- 
pletion altogether is also shown. It follows the Gaussian 
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FIG. 12: The final phase space distribution at the u = 30 
km/s sphere. In understanding this figure, it may help to take 
a look at the u = 30 km/s line of figure Q which corresponds 
to the central horizontal plane of this figure. The large red 
region to the right corresponds to unbound orbits. To the left 
(backwards 30 km/s) the phase space density is very low, as 
expected from the results of the solar depletion calculations. 
The leftmost part of the large red area corresponds to Jupiter 
crossing orbits, which are filled with the same density as the 
unbound orbits. 



down to about 2.5 km/s where it drops to zero. This is 
in perfect agreement with the results of Gould 0, and 
we can see this agreement as a test that our numerical 
routines are performing as they should. Our full numeri- 
cal routines without solar depletion is a numerical imple- 
mentation Gould's analytical arguments about diffusion 
in the solar system and our results should thus (as they 
do) agree in this case. We also show our raw numerical 
result, which is the outcome of our full simulation with 
solar depletion included. It is significantly lower than the 
Gaussian estimate in this low-velocity region, but not as 
low as the conservative (or ultra conservative) view. The 
general argument above that the time scales of solar de- 
pletion and diffusion are not too different and that some 
WIMPs should remain thus turns out to be valid. Hence, 
solar depletion kills some of the WIMPs at low velocities, 
but not as many as one could have feared. Also shown 
in the figure is our best estimate of the velocity distribu- 
tion, which is the same as our raw numerical result, but 
modified at low velocities (below 2.5 km/s) to include the 
effect of the eccentricity of the Earth's orbit, which will 
be explained now. 

The diffusion effects included so far does not provide 
means of filling the extremely slow (u < 2.5 km/s) or- 
bits. Such processes arise when the eccentricity of the 
Earth's orbit is taken into account. This could be done 
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FIG. 13: The radial velocity distribution of Earth crossing 
dark matter at the Earth. The curves labeled conservative 
and ultra conservative are the contributions from unbound, 
as well as unbound plus Jupiter crossing orbits respectively. 
The dash-dotted curve displays the result of ignoring solar 
depletion. The blue solid line represents our best estimate, 
including the effect of the eccentricity of the Earth's orbit. 
The thin line is the raw result from the numerical routines. 
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FIG. 14: The phase space density F(u) at low velocities. The 
upper curve is the Gaussian distribution. The thin solid curve 
is the outcome of our numerical simulations. The thick solid 
blue curve is our best estimate, where a population at low 
velocities (below 2.5 km/s) has been added due to the eccen- 
tricity of the Earth's orbit (see text). The dashed line, for 
comparison, shows the distribution in the conservative view 
where only unbound plus Jupiter-crossing orbits are included. 
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in a way similar to ordinary diffusion, but since u would 
no longer be fixed even in the one planet case, the block 
diagonal one planet time evolution operator would be 
polluted with new, off diagonal blocks, making the dif- 
fusion problem much more complicated. However, the 
eccentricity of the Earth's orbit will mean that the Earth 
diffuse slightly differently in the different parts of it's or- 
bit. This will cause a mixing of spheres of different u and 
thus cause an effective diffusion in the u direction. The 
size of this effect can be estimated using Eq. (2.10) of 
Gould's paper [2l|. Evaluation shows that for extremely 
slow particle orbits, the time scales can be as fast as one 
tenth of those of the ordinary diffusion, while in most 
other cases they are far slower. It is therefore quite rea- 
sonable to ignore these effects in our diffusion treatment 
at higher velocities. For the u < 5 km/s region the time 
scales of u-diffusion is comparable to the time scales of 
solar depletion, which makes it reasonable to assume that 
the phase space density is a slowly changing function with 
respect to u which means that the sharp cutoff at u = 2.5 
km/s is not physical. To estimate the phase space den- 
sity at these very low velocities, the mean density in the 
u £ [2.5, 5] km/s region is calculated and used as a mini- 
mum density in the whole u G [0, 5] km/s region. Another 
approach could have been to relocate the already existing 
mass to fill up the u < 5 km/s region evenly. However, 
this would underestimate the density in the u > 2.5 km/s 
region. Fig. [21 compares the raw result of the full nu- 
merical simulations, with this new best estimate and the 
Gaussian. The conservative view is also shown for refer- 
ence. 

We have now focused on the low velocity region of the 
velocity distribution. In Fig. ^| we show the full veloc- 
ity distribution for large velocities. At low velocities our 
results are confined by the free space Gaussian and the 
focused free space Gaussian of the conservative view. Fo- 
cused here means that the distribution as seen at the 
Earth is somewhat larger than in the halo due to the 
fact that particles are focused when they fall into the so- 
lar potential well. Thus, at typical galactic velocities, the 
Gaussian is somewhat lower than both our best estimate 
and the conservative view. 



VIII. CAPTURE AND ANNIHILATION RATES 

In the previous section, we have seen that our new es- 
timate of the WIMP velocity distribution is, especially at 
low velocities, considerably lower than earlier estimates 
based on the Gaussian approximation 5]. We will here 
investigate how this new velocity distribution affects first 
the capture rates of WIMPs in the Earth and secondly 
the annihilation rates of WIMPs in the center of the 
Earth. In this section, we will keep the discussion general 
and in section ITXl we will investigate the effects for the 
neutralino as a WIMP dark matter candidate. 




WIMP velocity u at the Earth (km/s) 



FIG. 15: The radial velocity distribution of Earth crossing 
dark matter at the Earth, linear scale. The curves are labeled 
as in figure lT^l Most of the velocity distribution is unchanged 
by the considerations in this report. The major difference 
between the Gaussian and the other distributions is that lat- 
ter have fewer slow particles, due to the effects of the solar 
potential. 



A. A new estimate of the capture rates... 

Given the velocity distribution derived in the previous 
section, we can now calculate the capture rate in the 
Earth with this velocity distribution. We will use the 
full expressions for the capture rate as derived by Gould 
in 1^1, but will also compare with the usual Gaussian 
approximation (as derived in 6]), as that is what most 
people use to calculate the capture rates. 

The calculation of the capture rates for an arbitrary ve- 
locity distribution is given in Q , we will here only briefly 
outline how the calculation is done. 

We divide the Earth into shells, where the capture from 
clement i in each shell (per unit shell volume) is given by 
H[Eq. (2.8)] 

where f(u) is the velocity distribution (normalized such 
that Jg 00 f(u) — n x where n x is the number density of 
WIMPs |3^|. The expression fi^to) is related to the 
probability that we scatter to orbits below the escape ve- 
locity, w is the velocity at the given shell and it is related 
to the velocity at infinity u and the escape velocity v by 
w = V 'u 2 + v 2 . The upper limit of integration is a priori 
set to u max — oo, but we will see below that due to kinc- 
matical reasons we can set it to a lower value (Eq. 
below). If we allow for a form factor suppression of the 
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FIG. 16: The capture rate of dark matter. This figure shows 
the rate at which dark matter particles are captured to the in- 
terior of the Earth, for a scattering cross section of a — 10~ 42 
cm 2 . The Gaussian-no solar depletion model gives the highest 
capture. The curves labeled ultra conservative and conserva- 
tive are the contributions from unbound, as well as unbound 
plus Jupiter crossing orbits respectively. For masses above 
150 GeV, our new capture estimate is considerably lower than 
that the Gaussian model. The peaks at low WIMP mass 
correspond to the masses of the included elements. A dark 
matter halo density of px=0.3 GeV/cm 3 is assumed. 
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we can evaluate wfl vi (w) and arrive at the expression 
i[Eq. (A6)] 
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with rrii the mass of element i. The Heaviside step func- 
tion plays the role of only including WIMPs that can 
scatter to a velocity lower then the escape velocity v. To 
simplify our calculations we can drop this step function 
in Eq. I|56|) and instead set the upper limit of integration 





Atomic 


Mass fraction 


Element 


number 


Core 


Mantle 


Oxygen, O 


16 


0.0 


0.440 


Silicon, Si 


28 


0.06 


0.210 


Magnesium, Mg 


24 


0.0 


0.228 


Iron, Fe 


56 


0.855 


0.0626 


Calcium, Ca 


40 


0.0 


0.0253 


Phosphor, P 


30 


0.002 


0.00009 


Sodium, Na 


23 


0.0 


0.0027 


Sulphur, S 


32 


0.019 


0.00025 


Nickel, Ni 


59 


0.052 


0.00196 


Aluminum, Al 


27 


0.0 


0.0235 


Chromium, Cr 


52 


0.009 


0.0026 



TABLE I: The composition of the Earth's core and mantle. 
The core mass fractions are from |22ll [Table 4] and the mantle 
mass fractions are from (22^1 [Tabic 2]. 
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We also need the scatte ring cross section on element i, 
which can be written as [12] [Eq. (9-25)] 



a A 2 (m x mi) 2 (m x + m p ) 2 
p % (m x + rrii) 2 (m x m p ) 2 



(59) 



where Ai is the atomic number of the element, m p is 
the proton mass and a p is the scattering cross section on 
protons. 

We now have what we need to calculate the capture 
rate. In Eq. (|53l) we integrate over the velocity for our 
chosen velocity distribution. We then integrate this equa- 
tion over the radius of the Earth and sum over all the 
different elements in the Earth, 



C 



^ dV 



(60) 



The capture rates depend on the mass and distribution 
of the elements in the Earth. The most important ele- 
ments are iron, silicon, magnesium and oxygen, of which 
iron is by far most important for WIMP masses over 100 
GeV. We use the Earth density profile as given in [2(j 
and for the element distribution within the Earth we use 
the values given in 22] [Table 2 for the mantle and Table 
4 for the core]. These values are listed in Table [I] 

Fig. 1161 shows the calculated capture rates, to be com- 
pared with that of a Gaussian distribution, with the 
Earth in free space. The Gaussian distribution is the 
one of equation in section llll Al This common model 
is equivalent to having the Earth taking the place of the 
Sun and removing the solar system (this is how weak cap- 
ture into the Sun is usually calculated). The scattering 
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FIG. 17: The ratio between the capture in various models 
and that of a Gaussian distribution in free space. The figure 
displays the quotient of the weak WIMP capture rates in the 
Earth in various models, and the capture given in the case of 
the commonly used Gaussian distribution. 



cross section between the nucleons and the WIMPs deter- 
mines the normalization only, and was taken to be 10~ 42 
cm 2 in Fig.^| We also show the resulting capture rates 
in the conservative and ultra conservative view, where 
the cutoffs at about 710 GeV and 410 GeV are clearly 
seen. These cutoff masses are higher than those in Gould 
and Alam as we have used the full integration over the 
Earth and not the average properties as in [8j (see section 
IIII Bl for a discussion of these cutoff masses) . 

It is of course interesting to compare the calculated 
capture rate with that given by the commonly used Gaus- 
sian distribution. This is done in Fig. 1171 where we divide 
by the capture rate in the Gaussian approximation. We 
clearly see that below 100 GeV, the different calculations 
agree to within about a factor of two. At higher masses 
the suppression is almost an order of magnitude, but not 
as bad as the feared conservative or ultra conservative 
views. 



where the first term is the WIMP capture, the second 
term is twice the annihilation rate Ta = \CaN 2 and the 
last term is WIMP evaporation. The evaporation term 
can be neglected for WIMPs heavier than about 5-10 
GeV and since we are not interested in these low-mass 
WIMPs we can safely drop the last term in Eq. Ij61|l . If 
we solve Eq. ()61|) for the annihilation rate Ta wc get 

Ya^-C^I, r=^= (62) 

where r is the time scale for capture and annihilation 
equilibrium to occur. In the Sun, equilibrium will for 
many WIMP models have occurred and the annihilation 
rate is at 'full strength', Ta — \C. In this case the 
annihilation rate is directly proportional to the capture 
rate. However, in the Earth, equilibrium has often not 
occurred, and we will have the more complex relation be- 
tween capture and annihilation rate, Eq. (|62|l . In the next 
section, we will show this for an explicit example, the neu- 
tralino in the Minimal Supersymmetric Standard Model 
(MSSM). Before looking at specific MSSM models, let's 
analyze Eq. (|62() to see the general trends. Let's denote 
the capture and annihilation rates in the usual Gaussian 
approximation by C G and T G respectively, whereas our 
new estimates are denoted C and Ta- Using the fact 
that the constant C a is the same in both scenarios, we 
can then write 

T A C ^ (Vgjj ^ [ & ;*0»r 
T A C ° tanh 2 (^) "l(^) 2 ;i «r 

(63) 

Hence, if equilibrium has occurred, the annihilation rate 
(and thus the neutrino-induced muon fluxes) are sup- 
pressed with the same factor as the capture rates, but 
if equilibrium has not occurred, the annihilation rate is 
suppressed with the square of the capture rate suppres- 
sion factor, i.e. the suppression is amplified. 

IX. APPLICATION TO THE 
SUPERSYMMETRIC NEUTRALINO 



B. ...and the annihilation rates 

We have seen that our new estimate of the capture rate 
in the Earth is, especially at higher masses, considerably 
lower than the usual estimate based on the Gaussian ap- 
proximation Since the neutrino-induced muon rates 
do not directly depend on the capture rate, but instead 
on the annihilation rate, we will here investigate how the 
annihilation rates are affected. 

The evolution equation for the number of WIMPs, N, 
in the Earth is given by 

^ = C - C A N 2 - C E N (61) 



So far, we have discussed the effects of our new esti- 
mate of the velocity distribution in general terms. We 
have seen that our estimate of the velocity distribution 
is significantly different from previous estimates at low 
velocities. We have also seen that the capture rates, es- 
pecially at higher WIMP masses are significantly reduced 
with a factor C/C . Hence, the annihilation rates (and 
the expected neutrino-induced muon fluxes) are reduced 
by a factor that lies in between (C/C G ) 2 and C/C G . We 
now want to investigate this suppression factor further 
and analyze the effects on the neutrino-induced muon 
fluxes. For this we need an explicit WIMP candidate. 
We will here assume that the WIMP is the lightest neu- 
tralino, that arises as a natural dark matter candidate 
in supersymmetric extensions of the standard model. In 
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the next subsection, we will briefly go through the su- 
persymmetric model we work in and will then continue 
to investigate the effects of our new velocity distribution 
on the annihilation rates and the neutrino-induced muon 
fluxes. 



~ 10' 



10 ' 



10 



24 



* °u > 10 km" 2 yr" 1 
4 oM < 10 km" 2 yr" 1 



J. Lundberg and J. Edsjd, 2003 



A. The neutralino as a dark matter candidate 

We will assume that the WIMP is the lightest neu- 
tralino in the Minimal Supersymmetric Standard Model 
(MSSM), i.e. the lightest neutralino, x°, is defined as the 
lightest mass eigenstate obtained from the superposition 
of four spin-1/2 fields, the Bino and Wino gauge fields, 
B and W 3 , and two neutral CP-even Higgsinos, H® and 
H° 2 : 
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NuB + N 12 W 3 + N 13 H° + N U H° 



(64) 



For a recent review of the MSSM and the neutralino as 
a dark matter candidate, see |2^|. The parameters of 
our phenomenologically inspired MSSM model are the 
Higgsino mass parameter ji 1 the gaugino mass parame- 
ter M2, the ratio of the Higgs vacuum expectation val- 
ues tan/3, the sfermion mass scale Mg, the mass of the 
CP-odd Higgs boson m,A, and the trilinear couplings for 
the third generation squarks A t and At,. We have made 
extensive scans of these parameters and have currently 
about a couple of hundred thousand models in our model 
database. 

For our actual calculations we use the DarkSUSY pack- 
age [24|. We only select those models that do not vio- 
late present accelerator bounds. The neutralino natu- 
rally has a relic density in the right ballpark, and we will 
further restrict this range by selecting only models with 
a relic density in the range 0.05 < VL x h 2 < 0.2. This 
range is a bit larger than the current best estimates 0, 
but to be conservative we choose work with this larger 
range. When calculating the relic density, we have in- 
cluded coannihilations between neutralinos and charginos 
(coannihilations also with sfermions in the MSSM is the 
subject of a future publication). 



B. Neutralino capture and annihilation 

We will here investigate how the annihilation rates are 
affected for specific MSSM models. In Fig. ED we show 
typical equilibrium time scales, r, for a set of MSSM 
models. As seen, the typical equilibrium time scales are 
much longer than the age of the solar system, t Q ~ 4.5 • 
10 9 years, and hence equilibrium has often not occurred 
in the Earth. 

As equilibrium has not occurred in the Earth, we can 
use Eq. I|62[l to see how the decrease in C will affect Fa- 

In Fig. El we show, for a set of MSSM models, how 
the annihilation rates are decreased. We also show the 
limiting cases for f Q 3> r and -C r. We can clearly see 



10 



22 



cr 



10 ' 



10 



10 



20 



19 



10 ' 



++ +++++++++++ 
+++ ++++++++ + 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 
++++++++++++++ 



+++++++++++++++++++++++++ 



10 



17 



'•"••••••"^ge of solar system 



10 ' 



10 



10 10 10 

Neutralino Mass (GeV) 



FIG. 18: The equilibrium time scales r for a set of MSSM 
models. As seen, all our models have equilibrium time scales 
longer than the age of the solar system. The green dots are 
models for which the neutrino-induced muon fluxes (with the 
usual Gaussian approximation) are larger than 10 km -2 yr -1 
and the blue plus signs indicate models with smaller fluxes. 
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FIG. 19: Fa/F^ versus the neutralino mass m x . The limiting 
cases for 4q >t and i© r are indicated in the figure. Most 
models have annihilation rate suppressions close to the lower 
curve since equilibrium has most often not occurred in the 
Earth. 
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that for most models, as equilibrium has not occurred, we 
are close to the (C/C ) 2 suppression of the annihilation 
rates. 



C. Neutrino-induced muon fluxes from the Earth 

So, given our calculated suppression of the annihila- 
tion rates, the neutrino-induced muon fluxes will also be 
suppressed by the same amount. We now ask ourselves if 
this suppression is too big to make the neutrino-induced 
muon fluxes too low to be observable in the MSSM. In 
Fig. 1201 we show in the left panel the neutrino-induced 
muon fluxes with the old Gaussian approximation. In the 
right panel, we show the neutrino-induced fluxes with our 
new estimate of the WIMP velocity distribution. We also 
indicate current limits from neutrino telescopes (Baksan 
153, Macro Amanda [2^| and Super-Kamiokande 

[29j) and anticipated sensitivities for future neutrino tele- 
scopes like IceCube [3(3] • Note that the IceCube limit 
shown here is a probably too optimistic, but we show it 
as limiting case beyond which a 1 km 3 neutrino telescope 
will not reach. For comparison, we also indicate the cur- 
rent direct detection limit by the Edelweiss experiment 
p^j . Models that are excluded by Edelweiss are indicated 
by green circles, whereas models that are not excluded 
are indicated with blue crosses. 

Comparing the left and the right figure, we clearly see 
that there is a significant suppression of the rates above 
about 100 GeV and above about 2000 GeV, the fluxes are 
too low to be observable even with future detectors. In 
the range between 100 GeV and 2000 GeV, where future 
neutrino telescopes still have a chance to detect a signal 
from the Earth, the prospects for doing so is clearly di- 
minished with our new estimate of the fluxes. Especially 
if one considers that all of the observable models in that 
range are already excluded by direct detection experi- 
ments. Note, however, that the comparison between di- 
rect detection and neutrino telescopes that we have done 
here is for a Maxwell-Boltzmann velocity distribution. 
As direct detection experiments are primarily sensitive 
to the high velocity tail of the distribution, whereas neu- 
trino telescopes are sensitive to the low velocity tail, the 
correlation between the two signals need not be as large 
as indicated in Fig. 1201 for a more realistic distribution. 
Below 100 GeV, the neutrino signal from the Earth is 
not reduced much with our new velocity distribution. In 
this range, neutrino telescopes are also in general more 
sensitive than direct detection experiments. 



X. CONCLUSIONS 

We have made a new estimate of the velocity distri- 
bution of WIMPs at the Earth due to diffusion in the 
solar system. We have included gravitational diffusion 
due to the Earth, Venus and Jupiter and depletion due 
to solar capture. Compared to the standard approxima- 



tion (i.e. that the solar diffusion can be approximated by 
the Earth being in free space and seeing the unperturbed 
Gaussian halo velocity distribution), our estimate is sig- 
nificantly lower at low velocities (below about 70 km/s). 
The main reason for this is that solar capture diminishes 
the WIMP population at these low velocities. If it were 
not for solar capture, our results would confirm the re- 
sults of Gould Q, i.e. that the velocity distribution as 
seen at the Earth is close to that we would see if the 
Earth was in free space. The diffusion effects of Jupiter, 
Earth and Venus would make the distribution look Gaus- 
sian, apart from a hole in the distribution below 2.5 km/s. 
This hole would however be filled due to the eccentricity 
of the Earth's orbit. However, solar capture suppresses 
the velocity distribution by about an order of magnitude 
at low velocities and this suppression propagates into a 
suppression of the same order of magnitude in the cap- 
ture rate. 

Since the annihilation rates depend on the capture 
rates, the annihilation rates are also suppressed. The 
amount of suppression, however, depends on if capture 
and annihilation is in equilibrium or not. If it is in equi- 
librium the annihilation rate suppression is the same as 
the capture rate suppression, but if we are far from equi- 
librium, the annihilation rate suppression is equal to the 
capture rate suppression squared. 

For one of the prime WIMP dark matter candidates, 
the neutralino in the Minimal Supersymmetric Standard 
Model (MSSM), it turns out that these are typically not 
in equilibrium and thus the annihilation rate suppression 
is equal to the capture rate suppression squared. The 
net result is that the annihilation rates will start being 
suppressed above about 100 GeV, and reaches a maximal 
suppression of about 10 -2 at around 1 TeV. Above about 
2 TeV, the expected fluxes are so low that future neutrino 
telescopes will not have enough sensitivity to see these. 

Finally, a word of caution should be applied to the in- 
terpretation of these new results. Even if we have done 
what we can to make sure that our new estimate is cor- 
rect, there are still approximations done and numerical 
uncertainties that need to be considered. E.g., in princi- 
ple one would like to do a full numerical simulation of the 
full diffusion process with an arbitrary halo distribution 
as input. That is not numerically feasible to do so in- 
stead we have relied on numerical simulations for the so- 
lar capture and on analytical calculations and arguments 
for the diffusion process. These analytical calculations 
are approximations with the aim to describe the diffu- 
sion processes correct in average. We think that these 
approximations are reasonable, but one should keep in 
mind that there are uncertainties involved in these ap- 
proximations. At higher masses, above about 1 TeV, we 
are very sensitive to the very details of the velocity dis- 
tribution at very low velocities (a few km/s). We have 
assumed that the eccentricity of the Earth's orbit fills the 
hole below 2.5 km/s. If this would not be the case, the 
suppression for high masses would be even larger than 
depicted here. 



21 



I 10' 

a 

a 

s 

e 

* 10 
10 

1 

10 
10 



J. Lundberg and J. Edsjd, 2003 



| 

c 
o 
s 



10 



E l " = 1 GeV 



0.05 < Oh^ < 0.2 



AMANDA 97 data, 02= 

BAKSAN 1997 

MACRO 2002 
- - SUPER-K 2002 prel : 

IceCube Best-Case I 



• •••• 

•••••• 

•••••• 

••••••••• 

••••••••• 



• • • »V» m ■ 

• •••••••••••• 

••••••••••••••••• 

••••••••••••••••a* 

• •••••••< ++ + 

•••••••••••••••• 

• •• ••••••••*- + + + + + + + 



- a CI > a. 



++ + 
++ + 

+++ +•••••••••••+••+++• 

+++ +++++++++++++++++++ 
++++++++++++++++++ ++++ 

-Kt-t + + + + + + + + + + + + + + + + + + 
i + i + i-F-TfTfTTTT-T-l-T-V ^1-"+"+"+"+" 

++++++++++++++++++ 
++++++++++++++++++ 

++++++++++++++++++ +++++ 
+ + + + + .+ + + + + + + + + ' ' 

++++++++ + + + + + + 

+ + + + + + + + + + + + + + + 
f+++++++++++++++ 
f +++++++++++ ++++ 
f++++++ ++++ +++++ 
f ++++++++++ ++++ 
f+++++++ +++++ ++++ 



+++++++ 
++++++ 



++++++ 
++++++ 
+++++++ 

+++++++ 



+++++ ++ 

+++++++++ 

+++++++++ 



+++++++++ 



++++++++ ++++++++++++++++++++++ 
+++++++++++++++++++++++ ++++++ 
+++++++++++++++++++++++++++ ++ 



++++++++++++++++++++++++++++ 
+l+-lt+l+H-H-l4 + + + + -lt + + -lf + -H + lH--lt-lt + ti- + + + 



10 



10 



10 



Neutralino Mass (GeV) 



9 

a 

c 

o 

3 



10' 



10 



J. Lundberg and J. Edsjd, 2003 



1 io' 



10 



10 



10 



10 



0.05 < Oh* < 0.2 



10 



AMANDA 97 data, 02 

BAKSAN 1997 

MACRO 2002 

SUPER-K 2002 prel : 

■ IceCube Best-Case 

New solar system diffusion 



••• ••••• 



•••••••• 



••••••• 

••••••• 

++++++ 



f >Oj 



+ +-H-+-k-L * « t f . 1 1 » ##. 
++++++•••••••••• 

++++++«♦♦••••••• 

++++++++++++++++ 
+++++ + ++++++++++ 

f++++++ ++++++++++++++ 
f+++++++++++++++++++++ 
f +++++++++++++ +++++•• 
++++++++++++++++++++•• 



SI 

+ ++++++ 
+++ +++++ 
++++++++++++++++++ +++ +++++ 



+++++++++++++++++ 

tltll tltltltltttttll tttt ll 



++++ ++++ 

I "tlill tfrttt t 



10 



10' 



10 



Neutralino Mass (GeV) 



FIG. 20: In the left panel we show the neutrino-induced muon fluxes in the standard Gaussian approximation, whereas in the 
right panel we show the fluxes based on our new estimate of the WIMP diffusion in the solar system. We also show the current 
limits of a few neutrino telescopes and an optimistic estimate for the future IceCube sensitivity. The current direct detection 
limit by the Edelweiss experiment j25| is also shown. Models that are excluded by Edelweiss are indicated by green circles, 
whereas models that are not excluded are indicated with blue crosses. 
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APPENDIX A: COMPARISON WITH 
FARINELLA'S CALCULATION OF NEAR 
EARTH ASTEROIDS 

This appendix considers the asteroids which fates were 
investigated by Farinella et al. 0. They calculated the 
fates of about 47 asteroids of which most are near Earth 
asteroids (NEAs). The result was that about a third of 
the considered asteroids were ejected to hyperbolic orbits 
or driven into the sun in less than 2 million years. This 
led Gould to consider the possibility that the population 
of gravitationally bound dark matter is heavily reduced 
by solar capture. This in turn lead to the Conservative 
and Ultra conservative views discussed in the introduc- 
tory sections. 

As a test of the Mercury integration package 01 > we 
have repeated the calculations of Farinella et al. using 



both Bulirsch-Stoer 16] and fifteenth-order Radao [l5|. 
These are quite complicated methods specially developed 
for solving the many-body problem. 

The actual fates of specific asteroids is of course de- 
pendent of the method used, and the accuracy param- 
eters of the calculation. Even with very high accuracy, 
convergence can not be expected since numerical errors 
propagate exponentially in chaotic systems. The initial 
conditions of our calculations are those of the online as- 
teroid database at U.S. Naval Observatory, epoch 11-22- 
2002. In addition to the time passed, some asteroids have 
been observed many times since 1994. However, one can 
still hope to imitate the general behavior by looking at a 
large set of initial values, regardless of what they repre- 
sent; asteroids or WIMPs. 

The results for the Bulirsh-Stoer method are presented 
in table [H] Of the 47 objects, four were ejected from 
the solar system and twelve were captured by the Sun in 
our calculations, whereas four were ejected and 14 were 
captured by the Sun in Farinella et al.'s calculations. We 
cannot expect to get exactly the same results on an indi- 
vidual basis, but are satisfied to see that we get roughly 
the same behavior as Farinella et al. 
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TABLE II: The asteroids integrated by Farinella et al. The 
numbers given are the times at which the asteroid collided 
with the Sun, or was ejected, in thousands of years. The (**) 
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